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ABSTRACT 

We present a new stellar dynamical mass measurement of the black hole in the nearby, SO galaxy NGC 3998. 
By combining laser guide star adaptive optics observations obtained with the OH-Suppressing Infrared Imaging 
Spectrograph on the Keck II telescope with long-slit spectroscopy from the Hubble Space Telescope and the 
Keck I telescope, we map out the stellar kinematics on both small spatial scales, well within the black hole 
sphere of influence, and on large scales. We find that the galaxy is rapidly rotating and exhibits a sharp central 
peak in the velocity dispersion. Using the kinematics and the stellar luminosity density derived from imaging 
observations, we construct three-integral, orbit-based, triaxial stellar dynamical models. We find the black hole 
has a mass of M BH = (8.1^°) x 10 8 M Q , with an /-band stellar mass-to-light ratio of M/L = 5.0^4 M Q /L Q 
(3(7 uncertainties), and that the intrinsic shape of the galaxy is very round, but oblate. With the work presented 
here, NGC 3998 is now one of a very small number of galaxies for which both stellar and gas dynamical 
modeling have been used to measure the mass of the black hole. The stellar dynamical mass is nearly a factor 
of four larger than the previous gas dynamical black hole mass measurement. Given that this cross-check has 
so far only been attempted on a few galaxies with mixed results, carrying out similar studies in other objects 
is essential for quantifying the magnitude and distribution of the cosmic scatter in the black hole mass - host 
galaxy relations. 

Subject headings: galaxies: active - galaxies: individual (NGC 3998) - galaxies: kinematics and dynamics - 
galaxies: nuclei 
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1. INTRODUCTION 

B lack holes reside at th e centers of most mass ive galax- 
ies (iMagorrian et al l 119981: iRichstone et"aTl 119981) and their 
masses (Mbh) are known to tightly correlate with properties 
of their h ost galaxy, like the bulge stellar velocity dispersion 
(ov: e.g.. iFerrar ese & M errittl 120001: iGebhardt et al.lE OQOb; 
Trema ine et al.l l2002t iGiiltekin et al.l l2009bt iGraham et alJ 
201 ll) and the bulge luminosity (Lhni; e.g., i Dresslerl 
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1989t iKormendv & Richstondll995t Marconi & Huntl 12003b 
Gtiltekin et al.ll2009bf) . Playing a fundamental role in estab- 
lishing the empirical Mbh _ host galaxy relationships over the 
past 15 years has been the Hubble Space Telescope (HST). In 
fact, HST spectroscopic observations have lead to about 50 of 
the 70 black hole mass measurements made to date. 

More recently, adaptive optics (AO) systems on large 
ground-based telescopes have become an attractive alterna- 
tive for obtaining high spatial resolution data. These systems 
correct for wavefront distortions caused by turbulence in the 
Earth's atmosphere using a natural guide star or a laser guide 
star (LGS) as a reference. They are often used in conjunction 
with an integral field unit (IFU), which can map out the two- 
dimensional (2D) velocity field very efficiently compared to 
traditional long-slit spectrographs. In addition, AO systems 
operating in the near infrared on 8 - 10m telescopes are able 
to probe deep into the central regions of galaxies that previ- 
ously presented significant observational challenges with the 
2.4m HST, such as giant ellipticals with low-surface bright- 
ness cores and galaxy nuclei obscured by dust. AO-assisted 
IFUs will undoubtedly have a major impact on the field of 
supermassive black hole detection, and black hole mass mea- 



surements using this technology have already become increas- 
ingly more prevalent in the literature (e.g.. iHoughtonetalJ 
| 2006t|Davies et alj|2006t iNeumaver et alj|2007t iNowak et all 
boidlRusli etal.1120111) . 

Measurements of black hole masses are most often made 
through the dynamical modeling of gas disks or stars. Gas 
dynamical modeling is conceptually simple, and the method 
can be applied to galaxies containing nuclear gas disks in cir- 
cular rotation. The main drawback of this technique is that 
the gas can be influenced by non-gravitational forces, and 
the assumption of circular rotation must be verified. Further- 
more, many of the observed gas disks in early-type galaxies 
exhibit velocity dispersions in excess of those expected from 
just rotational motion alone. This internal velocity dispersion 
can be quite large, ranging from ^100-500 km s" 1 (e.g.. 



van der Marel & van den Bosch 1998b; Verdoes Kleiin et all 
2002tlDalla Bonta et al.ll2009t) . The physical nature of the in- 
trinsic velocity dispersion is not understood. One interpreta- 
tion is that the intrinsic velocity dispersion is the result of lo- 
cal microturbulence, but the bulk motion of the gas remains a t 
the circular velocity (van der Marel & van den Boschll 998b). 
Alternatively, the intrinsic velocity dispersion could be dy- 
namically important, providing pressure support to the gas 
disk. Models that do not account for th is effect will un- 
derestimate the true black hole mass J^g.. lBarmet"aT]l200Tl 
INeumaver et al.ll2007tlWalsh et al.ll2010b . 

In contrast, stellar dynamical modeling is the most widely 
used method. However, extracting a secure black hole mass 
measurement from the observed kinematics can be challeng- 
ing due to the complexity of the models. Recent adjust- 
ments to several previous stellar dynamical measurements 
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have been made after identifying significant systematic ef- 
fects due to degeneracies between the dark matter h alo and 
stellar mass-to-light ratio dGebhardt & Thomas|[2009l) . insuf- 
ficient coverage of the p hase s pace occupied by tangentia l 
orbits (IShen & Gebhardtl 120101 ISchulze & Gebhardtl EoTTh 
and the effects of triaxiality (van den Bosch & de Zeeuw 
2010). Other possi ble sources of er ror, such as the choice 
of regularization ( Valluri et all 120041) and the uncertain in- 
clinations of e arly-type galaxies (e.g.. IVerolme et aflT 2002b; 
iGebhardt e t al. 2003; Krajnovic etal1 l2005l) . have lead to fur- 
ther debate about the accuracy of stellar dynamical mass mea- 
surements. 

Since gas and stellar dynamical modeling methods suffer 
from different systematic effects, carrying out consistency 
tests between the two techniques is crucial. Only by applying 
both mass measurement methods to the same object can we 
address important unanswered questions regarding the black 
hole scaling relations. Is there a systematic difference be- 
tween the masses derived from the two methods, and if so, 
how does this affect the slope of the relations? How much 
of the scatter in the Mbh - host galaxy relationships is the 
result of inconsistencies between stellar and gas dynamical 
measurements? Theoretical interpretations of Mbh _ host 
galaxy correlations range from the black h ole being an essen- 
tial component in ga laxy evolution (e.g.. ISilk & Reeslll998t 
iDi Matteo et alj 120051) to simply being the result of random 
mergers without the need for a coevolution of black holes and 
galaxies (lPengll2007t iJahnke & Macciol l201 ll) . Understand- 
ing the shape and cosmic scatter of the relations, especially at 
the sparsely populated upper- and lower-mass ends, are key to 
distinguishing between the various theories and for character- 
izing the black hole mass function. 

This necessary consistency check has only been at- 
tempted on a few galaxies with mixed results. Stel- 
lar an d gas dynamical models have been applied to IC 
1459 (VerdoesKleii netalJ l2000t ICappellari et alJ 12 002b) 
and NGC 3379 (IGebhardt et al.l l2000al: Ishapiro et al l 120061: 
Ivan den Bosch & de Zeeuwl l2010lh Unfortunately, in both 
these galaxies, the gas kinematics turned out to be disturbed 
and not compatible with regular disklike rotation. Conse- 
quently, the ionized gas does not give a useful measure- 
ment of the black hole mass, and the galaxies are not suit- 
able for a proper comparison between the two methods. In 
the cases of NGC 3227 and NGC 4151, the gas dynamica l 
black hole mass measurements by iHicks & Malkanl (120081) 
are gener ally consisten t with the s tellar dynamical ma sses re- 
porte d bvlDavies et al.1 (120061) and lOnken et al.1 (120071 ). How- 
ever, lOnken et al.l (|2007) could not find a single best fit- 
ting model for NGC 4 151, and label their s tellar dynami- 
cal result as tentative. Verd oes Kleiin et all d2002l) applied 
both methods to NGC 4335, and found a stellar dynami- 
cal mass that is at least five times larger than the gas dy- 
namical determination. The black holes in M87 and Cen A 
have been the subject of numerous gas and stellar dynami- 
cal studies. While early examinations of the stel lar dynam- 
ics (see review by [Kormendv & Richstone 1995, and refer- 
ences therein) produced similar mass estimates for the black 
hole in M87 as the gas dyn amical measurements (Har mset al.1 
ll994tlMacchetto et alJI 19971). t he m ost recent stellar dynami- 
cal mass from IGebhardt et al.l (1201 ll) is about a factor of two 
larger. Finally, there is a great variation in the Cen A gas 
dynamical black hole mass measurements, ma inly due to un- 
certai nties in the inclination of the gas dis k (Marco ni et al J 
1200 ll 2006, iHaring-Neu maver et al.l 120061, iKrainovic et alJ 



120071 INeumaver et al.l 120071). The la test gas dynamical mea- 
surement of INeum aver etall (120071) is si gnificantly smaller 
than the stellar dynamical measurement of ISilge et al.l (120051) 
(by about a factor of four), but is in agreement with th e stel- 
lar dynamical measurement by Cappellari et al. (2009). With 
only a few direct comparisons between the gas and stellar dy- 
namical techniques, clearly more objects need to be exam- 
ined. 

An excellent target for such a comparison study is NGC 
3998, which is a nearby, SO galaxy with a lo w ioniza- 
tion n uclear emission-line region (LINER) nucleus dHo et alJ 
119971) . The galaxy has a simple morphology suitable for 
both stellar and gas dynamical modeling. The nuclear 
gas kinematics have been previously studied using multiple 
slit positions of the HST Space Telescope Imaging Spec- 
trograph (STIS), and have been fit with a circ ularly rotat- 
ing thin-disk mod e l (Ide Francesco et al.1 12006). With the 
Ide Francesco et all (120061) measurement of 2.2 x 10 8 M Q 
(scaled to our assumed distance), the black hole sphere of in- 
fluence (r sp h ere = GMbh/c*) is expected to be 0."15, which can 
be resolved by IFUs combined with AO on 8- 10m ground- 
based telescopes. NGC 3998 also contains a bright, com- 
pact nucleus that is suitable for use as a LGS tip-tilt refer- 
ence. With th e bulge stellar velocity dispersion reported by 
iGultekinetaLl (l2009bl) of 305 km s~\ this galaxy is further 
interesting as falls at the upper end of the Mbh _ cr+ relation- 
ship. 

In this paper, we describe a measurement of the black 
hole in NGC 3998 using orbit-based stellar dynamical mod- 
els. We will compare the stellar dynamical black hole mass 
measurement to th e existi ng gas dynamical measurement by 
Ide Francesco et all {2006) in order to test whether the two 
methods give consistent results when applied to the same 
galaxy. We begin by presenting the imaging observations and 
the determination of the galaxy's stellar mass distribution in 
$2] and $3] We describe the spectroscopic observations in §|4] 
and the measurement of the stellar kinematics in §0 Mod- 
els of the point-spread function (PSF) are discussed in §|6] In 
$7] we provide a brief overview of the stellar dynamical tech- 
nique, and its application to NGC 3998, as well as present the 
modeling results in §|8] Finally, in §|9] we examine the stellar 
orbital structure of the galaxy, compare the stellar dynamical 
mass measurement to gas dynamical determination, and place 
the black hole in NGC 3998 on the M B h - host galaxy rela- 
tionships. Througho ut this paper, we adopt a distance to NGC 
3998 of 13.7 Mpc dTonrv et al.ll200ll iMei et al1l2007l) . 

2. IMAGING OBSERVATIONS 

Imaging observations are essential components in the orbit- 
based stellar dynamical models, and are used to derive the 
galaxy's luminosity density. We obtained a Wide Field Plan- 
etary Camera 2 (WFPC2) F791W image of NGC 3998 cen- 
tered on the Planetary Camera (PC) detector from the HST 
archive. The HST image of NGC 3998 was originally ac- 
quired under program GO-5924, and the total exposure time 
was 100 s. The WFPC2 image, with a pixel scale of 0"046 
and 0"1 for the PC and Wide-Field chips, provides high an- 
gular resolution imaging suitable for constructing a luminous 
mass model near the black hole, however the image only ex- 
tends out to a radius of ^10" (0.6 kpc). We therefore also 
included in the analysis a deep Ks-band image (courtesy of 
Lasker et al., in preparatio n) taken with the Wide-field In- 
fraRed Camera (WIRCam: IPuget et al.1 [20041) on the 3.6m 
Canada-France-Hawaii Telescope (CFHT). The image has a 
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spatial scale of 0"3 pixel" 1 and covers a 30' x 30' field. The to- 
tal exposure time was 442 s, and we were able to measure the 
galaxy's surface brightness profile out to a radius of ^250" 
(16.6 kpc). The ground-based image was used to generate 
the mass model on large spatial scales, which is helpful in 
constraining the intrinsic shape of the galaxy in the stellar dy- 
namical models. 

3. STELLAR MASS PROFILE 

We parameterized the WFPC2 F791W and WIRCam 
Ks-band images of NGC 3998 as the sum of 2D Gaus- 
sians using the M ulti-Gaussian Expansion (MGE) method 
(Cappellari 2002a). The MGE method has often been used 
to model high-resolu ti on and ground-based i mages (e.g., 
Verolme et all I2002ri JC appellari et all 120061: lOnken et al.1 
2007t lDalla Bo nta et alJl2009t iKrajnoyic et al.ll2009h . and we 
used the MGE method here because it has the advantage of 
being able to reproduce a large range of density distributions 
while also allowing for the deprojection to be carried out an- 
alytically. MGE models were fit to both images simultane- 
ously while also accounting for the HST PSF. The PSF itself 
was described as the sum of 25 positive and negative Gaus- 
sians, which were f ound by applying th e MGE software to 
a Tiny Tim model dKrist & Hookl 12004). Our best-fit MGE 
parametrization of the galaxy was composed of 12 Gaussians, 
where the innermost Gaussian was constrained to be round 
and the position angles (PA) of the Gaussian components were 
required to be the same. This model produced an excellent 
fit to the imaging data, as can be seen in Figure Q] In addi- 
tion, in Figure |2] we present the surface brightness profile, 
the MGE model, and the percent error averaged along the az- 
imuthal axis. The MGE model was corr ected for galactic ex- 
tinction using the Schleg eTet al.l d 19981) values given by the 
NASA/IPAC Extragalactic Database and the surface density 
was co nverted to /-band solar units using the WFPC2 calibra- 
tion by Dolphin!; (120001) . In Table [T] we provide the best-fit 
values of the MGE parameters. We further note that there 
are no significa nt emission lines within the F791W bandpass 
dHoetalJll993l) . 

Although our adopted MGE model requires each Gaussian 
component to have the same PA, we fit an additional MGE 
model that allowed for isophotal twists. We detected very 
small changes to the PA, of typically < 1°, between the com- 
ponents. Furthermore, when using this MGE parameteriza- 
tion as input into the orbit-based stellar dynamical models, we 
found worse agreement between the observed and predicted 
stellar kinematics. As a result, we do not consider the MGE 
model that allows for isophotal twists any longer, and focus 
solely on the MGE parametrization that constrains each com- 
ponent to have the same PA. 

When constructing mass models of galaxies containing an 
active galactic nucleus (AGN), often the innermost Gaussian 
of the MGE model is assumed to arise from non-thermal 
emission and is excluded from the stellar mass distribution. 
The nucleus of NGC 3998 h as been spectro scopically clas- 
sified as a Type 1.9 LINER dHo et alj|1997h . and an unre- 
solved, variable UV source, a compact radio source, and an 
X-ray source have all been detected at the galaxy's center 
dHummel et al.lll984tlFabbiano et al.lll994lMaoz et alJl2005h 
iRoberts & Warwickll2000t iPellegrini et al.ll2000h . All of this 
evidence suggests that NGC 3998 hosts an AGN. However, 
the galaxy also exhibits a cuspy surface brightness profile, 
and some starlight may still be contained within the inner- 
most Gaussian. We therefore ran orbit-based stellar dynami- 
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FIG. 1. — Isophotes of the MGE model (black) are compared to the inner 
regions of the HST WFPC2/PC F791W image (top) and the CFHT WIRCam 
Ks-band image (bottom) of NGC 3998. Contours are logarithmically spaced, 
but arbitrary. 

cal models using MGE parameterizations that both included 
and excluded thi s cen tral Gaussian, and will discuss the re- 
sults in $8] an d 3l8.il These descriptions represent the two 
extremes: in one case all the light from the innermost Gaus- 
sian is attributed to the stars and in the other scenario all the 
light is assigned to the AGN. 

4. SPECTROSCOPIC OBSERVATIONS 

We observed NGC 3998 using the IF U OH-Suppressing 
Infrared Imaging Spectrograph (OS IRIS: lLarkin et al.ll2006h 
assisted by the LGS AO system (Wizinowich et al. 2006; 
Ivan Dam et alj f2006) on the 10-m Keck II telescope. The 
OSIRIS data provided high angular resolution observations 
of the stellar kinematics over a 2D field and allowed us to 
resolve r sp h er e- We also retrieved archival STIS observations 
of NGC 3998 from the HST archive. The STIS observations 
have a well-characterized PSF and provided us with additional 
high-resolution data of the stellar kinematics along a single 
PA. Preliminary work from a stellar dyn amical analys i s usin g 
this same STIS data was described by iBower et al.l (|2000). 
Past studies have shown the importance of acquiring large- 
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FIG. 2. — Top panel: Comparison between the observed surface brightness 
profile (open squares) and the best-fitting MGE model (solid line) averaged 
along the azimuthal axis. Individual components making up the MGE model, 
which have been convolved with the HST PSF, are also shown. Bottom panel: 
The percent error between the model and data. 



TABLE 1 
Best-fit MGE Parameters 



j 

(1) 


log Ij (L Q ,/ pc 2 ) 
(2) 


«rj (") 
(3) 


«5 
(4) 


1 


6.445 


0.017 


0.990 


2 


5.405 


0.078 


0.940 


3 


4.808 


0.223 


0.962 


4 


4.445 


0.504 


0.946 


5 


3.814 


1.037 


0.756 


6 


4.096 


1.200 


0.903 


7 


3.973 


2.042 


0.881 


8 


3.641 


4.134 


0.799 


9 


3.347 


7.223 


0.813 


10 


2.444 


21.210 


0.796 


11 


2.149 


34.812 


0.776 


12 


1.034 


73.492 


0.866 



NOTE. — Column (1) gives the Gaussian com- 
ponent number, column (2) is central surface bright- 
ness of each Gaussian, column (3) provides the dis- 
persion of Gaussian along the major axis, and col- 
umn (4) lists the axis ratio. All of the Gaussian 
components were constrained to have the same posi- 
tion angle, which was determined to be -41.5°. The 
primed variables indicate projected quantities. Stel- 
lar dynamical models were nan using MGE models 
that both include and exclude the central component 

0' = D- 



scale observations of the kinematics in order to properly con- 
strain t he orbital distribution in the stellar dy namical models 
(e.g., [Shapiro et al. 2006; Krai novic et alj|2009h . For this rea- 
son, we also obtained long-slit observati ons using the Lo w- 
Resolution Imaging Spectrograph (LRIS: IOke et al Jfl 9951) on 
the 10-m Keck I telescope. In the sections that follow we 
describe the OSIRIS, STIS, and LRIS observations and data 
reduction. 



4.1. LGS AO OSIRIS Observations 

The OSIRIS observations were acquired on 2009 May 1 
and 2 during the first half of each night. We used the 0."05 
spatial scale and the broadband K (Kbb) filter, providing cov- 
erage of 1.965-2.381 /im. We aligned the long axis of the 
IFU with th e major axis of the nucle ar emission-line disk at 
PA = 308° (Ide Francesco et al.ll2006h . Using the nucleus of 
NGC 3998 as the tip-tilt reference, we were able to achieve 
a good AO correction and we estimate that the PSF core 
had a full width at maximum intensity (FWHM) of ~0."1 
(see $6]l. We obtained 300 s exposures of the galaxy nu- 
cleus and sky, following the sequence Object-Object-Object- 
Sky Object-Object-Object-Sky, and recorded 3.8 hours of on- 
source integration time. Each object exposure was dithered 
0."2 perpendicular to the long-axis of the IFU in order to re- 
move bad pixels and to obtain a slightly wider field of view 
of 1"2 x 3 "2. In addition, we observed K and M giant stars 
for use as velocity templates in the same observational setup 
as was used for the galaxy, and AO V stars for telluric correc- 
tion. 

The data were reduced using the OSIRIS data reduction 
pipeline v2.33 The pipeline includes subtraction of a sky 
frame, glitch identification and cosmic ray cleaning, extrac- 
tion of the spectra into a data cube (with two spatial dimen- 
sions, x and y corresponding to the short and long axis of the 
IFU, and one spectral dimension), wavelength calibration, at- 
mospheric dispersion correction, and telluric correction us- 
ing the extracted one-dimensional (ID) spectrum of an AO V 
star. While the reduction pipeline is capable of performing 
scaled sky subtraction following the algorithm outlined by 
iDaviesI d2007l) . we found that the current implementation did 
not work well for our data set and produced strong sky resid- 
uals in the reduced cubes. We instead opted for a direct sub- 
traction method, where we subtracted a sky exposure from 
the previous two object frames and from the following object 
frame. Thus, the sky and object frames were separated in time 
by at most 600 s, and the direct subtraction method worked 
well. After all of the individual galaxy cubes were reduced, 
we determined the spatial x and y offsets between the first 
galaxy cube and each of the remaining object cubes. This was 
done by summing along the wavelength axis of the data cubes 
to generate flux maps, then cross-correlating the images. Us- 
ing these offsets, the individual object cubes were aligned and 
averaged to produce the final galaxy data cube. The velocity 
template star observations were reduced in a similar manner, 
but an additional step was included to extract a ID spectrum 
from the data cube using a circular aperture with a radius of 7 
pixels. 

4.2. Archival STIS Observations 

We retrieved STIS observations of NGC 3998 made on 
1999 March 10 under program GO-7350 from the HST 
archive. NGC 3998 was observed with the 52x0.2 aper- 
ture at eight slit positions, all of which were aligned at a PA 
of 152°, within 15° of the galaxy's major axis. Two succes- 
sive exposures were obtained at each slit position to aid in the 
removal of cosmic rays, and a shift of about four pixels in the 
spatial direction was made between positions. The average 
exposure time recorded at each slit position was 2730 s. The 
G750M grating, centered on 8561 A, was read out in an un- 
binned mode, which provided a scale of 0."0507 pixel" 1 along 

5 http://irlab.astro.ucla.edu/osiris/pipeline.html 
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the spatial axis and 0."554 A pixel" 1 along the dispersion di- 
rection. Following the observation of the galaxy at each slit 
position, a contemporaneous fiat field was taken. 

In addition to the galaxy observations, the K3 III star HR 
260 was observed as part of Program 7350 for use as a ve- 
locity template. In order to build up a larger library, we also 
searched the HST archives for other suitable stars, using as 
a guide the previously published stellar dynamical black hole 
mass measurements made from STIS data. We retrieved an 
additional four stars, observed using the same STIS setup, 
from the archive: the K0 III star HR 7615 (GO-7566), the 
G8 III star HD 141680 (GO-8591), and the K2 III and G5 V 
stars HD 73471 and HD 115617 (GO-8928). 

The data were reduced by running individual IRAF0 tasks 
within the standard Space Telescope Science Institute (STScI) 
CALSTIS pipeline. With these tasks, the overscan region was 
trimmed, the bias and dark files were subtracted, flat-field 
corrections were applied, CR-split exposures were combined 
to reject cosmic rays, wavelength and flux calibration were 
performed, and the data were rectified for geometric distor- 
tions. We chose to run individual segments of the pipeline 
independently, followi ng the procedure detailed in the STIS 
Data Handbook v6.0 dBostroem & Prof fittl 1201 1). instead of 
running the wrapper CALSTIS program so that a couple of 
important modifications could be made. 

Specifically, we found that the dark calibration file provided 
by STScI, which is constructed from several long dark expo- 
sures taken over the course of a week, contained very strong 
hot pixels that did not subtract well and left numerous neg- 
ative pixels in the galaxy spectral images. This poor sub- 
traction is likely due to pixels that vary on timescales of less 
than a week . We therefore followed a similar approach to 
iBower et al.l (1200 ll) . and modified the dark file by replacing 
any pixels that deviated by more than 8rr with a median value. 
After the initial processing was complete, we in cluded an ad- 
dition al cleaning step and used LA-COSMIC ( Ivan Dokkuml 
12001b to remove the hot pixels remaining after the subtraction 
our modified dark file. We then applied the geometric distor- 
tion correction to the cleaned images. Additionally, we used 
the contemporaneous flat field exposures taken at each of the 
slit positions to remove the obvious fringe pattern affecting 
the spectra. IRAF tasks were used to reduce and normalize 
the flat field observations, shift and scale the fringes in the 
normalized flat such that they would match those in the galaxy 
images, and to finally apply the fringe correction. 

The reduced, geometrically rectified, fringe-corrected spec- 
tral images at each of the eight slit positions taken along the 
galaxy's major axis were then aligned and combined to pro- 
duce the final 2D galaxy spectrum. The same procedure was 
applied to the five velocity standard stars observed with STIS, 
but as a final step, we constructed a ID spectrum by adding 
together 3 spectral rows above and below the central row from 
the 2D images. 

4.3. Long-Slit LRIS Observations 

The LRIS observations were obtained during the first half 
night on 2009 April 15. On the red side, we used the 831 lines 
mm -1 grating centered on 8200 A with a l"-wide slit, produc- 
ing a scale of 0.92 A pixel" 1 in the dispersion direction and 
0"21 1 pixel" 1 in the spatial direction. We placed the long-slit 

6 IRAF is distributed by the National Optical Astronomy Observatory, 
which is operated by the Association of Universities for Research in Astron- 
omy under cooperative agreement with the National Science Foundation 



along four P As: along the major axis of the nuclear gas disk 
(PA = 308°; de Francesco et al. 2006), along the minor axis 
(PA = 218°), and at two intermediate angles (PA = 353° and 
PA = 263°). At each PA we acquired a 2 x 600 s observation, 
with the exception of the minor axis PA, where we recorded 
a 2 x 300 s observation. Between individual exposures, we 
dithered by 90" along the length of the slit. We additionally 
observed K giant stars for use as velocity templates, and the 
flux standard star Feige 34, in the same observational setup as 
was employed for NGC 3998. 

We reduced the LRIS data using IRAF. The initial process- 
ing steps of the 2D spectral images included trim and bias 
subtraction, flat fielding, and cosmic ray cleaning. For the 
NGC 3998 exposures, we also geometrically rectified the 2D 
spectra so that both the wavelength and spatial axes would 
run parallel along rows and columns, respectively. During 
this process, we used the Hg, Ar, and Ne arc lamp exposures 
acquired immediately following the NGC 3998 observations 
to wavelength calibrate the spectral images. We removed the 
sky background by subtracting the two dithered spectral im- 
ages taken at each PA from one another. The difference im- 
ages were spatially aligned, averaged together, and trimmed 
to remove the negative features resulting from the pair sub- 
traction. As a final reduction step, we compared a calibrated 
ID standard star spectrum to a ID spectrum of the galaxy ex- 
tracted using the APALL task to determine the multiplicative 
factor (as a function of wavelength) that was needed to bring 
the two spectra into agreement. The scaling was then applied 
to the 2D galaxy spectrum to produce a final, flux-calibrated 
spectral image of NGC 3998 at each PA. 

The velocity template stars were reduced following a 
slightly different procedure because we only required a ID 
spectrum from the 2D image. After the initial image process- 
ing, we applied the same spectral rectification and wavelength 
calibration to the 2D spectra of the velocity template stars as 
was used for the NGC 3998 spectra. We performed a sec- 
ondary correction to the wavelength calibration by applying a 
small linear shift along the dispersion axis so that two bright 
sky lines bounding the Ca II triplet region matched the known 
wavelengths of 8430 A and 8886 A. We then extracted a ID 
spectrum with APALL using a 5 "3 aperture radius and sky re- 
gions that were 8 "4 in width beginning at a distance of 23". 
The velocity template spectra were flux calibrated using an 
extracted ID spectrum of the standard star Feige 34. 

5. MEASURING THE STELLAR KINEMATICS 

From the OSIRIS, STIS, and LRIS data, we extracted 
ID spectra over a range of spatial locations and measured 
the line-of-sight velocity distribution (LOSVD), which we 
parameterized with the first four Gauss-Hermite moments 
(V, a, hj, and /i4). Measurement of the LOSVD requires 
spectra with a high signal-to-noise ratio (S/N). Typically, 
a S/N of ^40 per spectral and spatial resolution element 
are needed in order to reliably determine I13 and /14, which 
quantify as ymmetric and symmetric dev i ations from a Gaus- 
sian (e.g.. Ivan derMarel & Franxl 119931 iBenderet all 119941 
IStatler I ft99$. iKrainovic et alj 120091). We used the Voronoi 
binning algorithm (Cappellari & Copin] 120031) in order to 
achieve the best spatial resolution possible given such a re- 
quirement on the minimum S/N. The stellar kinematics in 
each of the spatial bins were then m easured using the penal- 
ized p ixel fitting (pPXF) method of Cappellari & Emsellem 
(2004). With this method, logarithmically rebinned galaxy 
spectra are fit in pixel space using a stellar template that is 
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convolved with the LOSVD. 

The pPXF code minimizes template mismatch by finding an 
optimal stellar template that is composed of a linear combina- 
tion of spectra from a library. The libraries were designed to 
contain stars that are representative of the galaxy's expected 
stellar population. When measuring the stellar kinematics 
from the OSIRIS data, we used a template library composed 
of 20 stars, consisting of mostly K and M giants (spectral 
types K0 - K5 and MO - M5), a K4 V star, and a G9 V star. 
The STIS kinematics were determined using a library of five 
stars (G8 III - K3 III, and G5 V), and the LRIS kinematics 
were extracted using a library of eight K0 III - K3 III stars. 

The optimal stellar template was determined just once by 
using pPXF to fit a very high S/N spectrum generated by sum- 
ming the spectra over all spatial locations. Once the optimal 
stellar template was found, pPXF was used to measure the 
kinematics in each spatial bin by holding fixed the relative 
weights of the stars that comprise the optimal template. We 
did, however, allow the coefficients of an additive Legendre 
polynomial of degree 1 and a multiplicative Legendre poly- 
nomial of degree 2 to vary. The additive polynomial was used 
to model the AGN continuum, and the multiplicative poly- 
nomial was included to correct continuum shape differences 
between the optimal stellar template and the galaxy spectra 
due to imperfect flux calibration and reddening. 

Errors on the kinematic measurements were estimated us- 
ing a set of Monte Carlo simulations. During each realization, 
random Gaussian noise was added to the observed spectrum 
from each bin based upon the standard deviation of the pPXF 
model residuals. The LOSVD was measured with the penal- 
ization turned off in order to produce realistic uncertainties. 
From 100 realizations, we determined the la errors from the 
standard deviation of the resulting distributions. 

We also tested the robustness of our kinematic measure- 
ments by assuming different continuum models, each with 
various combinations of additive and multiplicative polyno- 
mials of degree 1-3. We additionally used pPXF to measure 
just V and a in each spatial bin, as well as allowed for the rel- 
ative contributions of the stars making up the optimal stellar 
template to vary between spatial bins. Regardless of the con- 
tinuum model used, the number of Gauss-Hermite moments 
extracted, or whether a new optimal template was fit for each 
bin, we found consistent results for all but a few spatial bins. 
The kinematics measured from these few spatial bins were 
deemed unreliable and excluded from further analysis. 

As a final step, we determined the systematic offsets in 
the odd Gauss-Hermite moments directly from the data us- 
ing th e point-symmetrization rou t ine de scribed in Appendix 
A of Ivan den Bosch & de Zeeuwl d2010h . The offset for the 
first odd Gauss-Hermite moment corresponds to the recession 
velocity of the galaxy, and a systematic shift in the second 
odd Gauss-Hermite moment may result from slight template 
mismatch effects. The offsets were subtracted from the V and 
hj, values determined with pPXF. 

Although the systematic offsets in the odd moments 
were measured with the point-symmetrization code, we did 
not symmetrize the kinematics before using them as in- 
puts into the orbit-based dynamical mode ls. Often the 
observed kinematics are sym metrized (e.g., Gebhar dt et al.l 
2003; CappellarietaLl E006) bee ause the dynamical models 
can only produce predicted kinematics that are bi- or point- 
symmetric. Thus, symmetrizing the input kinematics reduces 
the noise in the observations and aids in visual comparisons 
between various models. However, the symmetrization rou- 
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FIG. 3. — The best-fit broadened optimal stellar template (red) is compared 
to the observed NGC 3998 spectrum from a single OSIRIS lenslet located 
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centered on x = —Of 1, y = — 1"0 (middle). Green dots show the fit residuals 
and have been shifted by an arbitrary amount. The optimal template used 
to measure the kinematics is dominated by a K5 III star, whose spectrum is 
shown in the bottom panel. 



tine assumes that the galaxy nucleus is exactly centered on 
one of the IFU lenslets. Whi le this is true for t he reduced data 
sets from the S AURON IFU dBacon et alJ2001l) . for which the 
symmetrization code was originally designed, the assumption 
does not hold for the OSIRIS data. Consequently, our best-fit 
stellar dynamical model, presented in §|8] was determined us- 
ing non-symmetrized kinematics, but in §!8.1l we also test the 
effect on Mbh of symmetrizing the kinematics before running 
the orbit-based stellar dynamical models. 

5.1. OSIRIS Kinematics 

We used the pPXF method to measure the kinematics from 
the OSIRIS data over a wavelength range of 2.22-2.38 fim. 
For NGC 3998, this wavelength region included three very 
strong CO bandhead absorption features: (2-0) I2 CO, (2- 
l) 12 CO, and (4-2) 12 CO. The kinematics were constrained 
mainly by these prominent CO bandheads, although weaker 
Ca I and Mg I absorption lines were also detected. We note 
that no emission lines were detected with the OSIRIS Kbb 
filter. The stellar kinematics were measured from 90 spatial 
bins, 10 of which contain just a single lenslet and fall within 
the central ~0."15 region. The S/N of the spectra (taken to 
be the ratio of median value of the binned spectrum to the 
standard deviation of the pPXF model residuals) ranged from 
39-73, with a median value of 58. The optimal template was 
composed of eight stars, with a K5 III star contributing a ma- 
jority of the flux. In Figure|3] we show example spectra from 
two spatial bins and the optimal template broadened by the 
best-fit LOSVD, as well as the spectrum from the K5 III star. 

The kinematics show that the galaxy is rapidly rotating, 
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with velocities of ±115 km s" 1 over the OSIRIS IFU, and 
there is an anti-correlation between hi and V. Although the 
hi, kinematics are expected to be symmetric about the center, 
with a peak or dip at the nucleus, we find h\ to be noisy, with- 
out any clear trends. In contrast, the velocity dispersion shows 
a very strong peak at the nucleus, rising from a ^270 km s" 1 
in the outer regions to a ^560 km s -1 at the center. The errors 
on the velocity dispersion measuremen ts ar e quite large at nu- 
cleus and will be discussed further in $8.11 but outside of the 
innermost lenslets, we find typical errors of 10 km s , 13 km 
s , 0.02, and 0.03, for V, a, hi, and h\, respectively. 

NGC 3998 contains an AGN, which dilutes the observed 
stellar features, and can lead to difficulties in extracting the 
kinematics near the nucleus. During the spectral fitting with 
pPXF, we included an additive polynomial to account for the 
AGN dilution. However, we measured large uncertainties on 
the kinematics extracted from four of the innermost lenslets, 
suggesting that pPXF cannot disentangle the AGN from the 
stellar component in a consistent manner. Our final best-fit 
model presented in §|8] was found using the kinematic mea- 
surements and errors from all the OSIRIS bins, however we 
further tested the effect on Mbh whe n the central kinematics 
are excluded from the modeling (see $8. It . 

5.2. STIS Kinematics 

We measured the stellar kinematics from the Ca II triplet 
lines using the pPXF routine to fit an optimal stellar template 
over the wavelength range 8420 - 8830 A. The optimal tem- 
plate was composed of two stars: a K2 III star, which con- 
tributed most of the flux to the template, and a K3 III star. 
Near the nucleus, there was substantial AGN contamination, 
which diluted the the Ca II triplet lines, as well as a clear de- 
tection of an emission line, which is most likely the [Fe II] line 
at a rest wavelength of 8618 A. During the spectral fitting, we 
excluded a 50 A wavelength region surrounding this emission 
feature to prevent biases when measuring the LOSVD. Fur- 
thermore, we were unable to measure the stellar kinematics 
within ^0."2 of the nucleus due to the AGN. In the end, we 
extracted reliable kinematics from eight spatial bins extend- 
ing from about 0''2 to 1" from the nucleus. The spectra in the 
eight bins had a S/N of 3 1 -40 with a median value of 38. We 
present the nuclear spectrum of NGC 3998, along with exam- 
ples of the spectral fitting for two bins from which the stellar 
kinematics could be measured in Figure |4] 

Like the OSIRIS kinematics, the STIS kinematics show that 
the galaxy is rotating quickly with velocities of ±100 km s" 1 
over the inner ±1", that there is a rise in velocity dispersion 
from a ~ 290 km s" 1 at 1" to a ~ 370 km s" 1 at 0."2, and that 
/13 is anti-correlated with the velocity. Again, we are unable to 
detect any significant trends in h\. The median errors on the 
STIS kinematics were 21 km s" 1 , 25 km s" 1 , 0.05, and 0.06 
for V, a, hi, and hi,, and are larger than the uncertainties on 
the OSIRIS kinematics due to lower S/N spectra. Given these 
uncertainties, the STIS data should not have a substantial im- 
pact on the Mbh determination, and we found that this was 
indeed the case when running tests that excluded the STIS 
measurements. However, for completeness, the final best-fit 
model described in ^8] was determined using the STIS data, 
along with the OSIRIS and LRIS measurements. 

5.3. LRIS Kinematics 

We used the pPXF method to measure the kinematics from 
the Ca II triplet lines over a wavelength range of 8480 A - 
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the spectral fitting. The bottom panel shows the STIS spectrum of a K2 III 
star, which contributed most of the flux to the optimal stellar template. 



8830 A. The kinematics in each spatial bin were determined 
using an optimal stellar template constructed from the linear 
combination of a K3 III and a K0 III star. Ultimately, the kine- 
matics were measured from a total of 170 spatial bins over the 
four PAs. The S/N in the bins ranged from 60-107, with a 
median value of 82. We present example fits to the galaxy 
spectra in addition to the K3 III template in Figure [5] 

We were able to measure the kinematics out to ~15" (^1 
kpc) along the major axis and along the intermediate angles 
of PA = 353° and PA = 263°. The kinematics along the minor 
axis were measured out to ^8". The large-scale kinematics 
exhibit features similar to those seen from the high-resolution 
OSIRIS maps and STIS data. The stars show rotation with 
V between ±200 km s" 1 , the velocity dispersion increases 
steeply toward the center to values of 350-400 km s" 1 , and hi 
is anti-correlated with V across the inner ~10". In addition, 
we find that h\ is symmetric about the center with a slight 
peak to values of 0.05 at the nucleus. The median errors for 
V, a, /13, and h\ were 8 km s" 1 , 9 km s" 1 , 0.02, and 0.03 for 
V, a, hi, and hi,. 
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shown in red, and the fit residuals are displayed in green. For comparison, 
we plot the LRIS spectrum of the K3 III star that contributes most of the flux 
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We find excellent agreement between kinematics measured 
from the OSIRIS, STIS, and LRIS data. While template mis- 
match is a common source of uncertainty when extracting 
stellar kinematics, at least it is encouraging to find consistent 
measurements from different instruments over different wave- 
length regions using different velocity template stars. In Fig- 
ure [6] we show the four Gauss-Hermite moments measured 
from the three different instruments as a function of projected 
distance from the nucleus. 

6. PSF MEASUREMENTS 

An important input into the stellar dynamical models is the 
spatial PSF. However, determining the OSIRIS PSF is espe- 
cially complicated because the quality of the AO correction 
varies with time and numerous data cubes (obtai ned over the 
course of two nights) were mosaiced together. iDavies et al.l 
(2004) highlight several ways to estimate the PSF from AO- 
corrected IFU data, and for this work we used a method that 
compares the OSIRIS data cub e to a deconvolve d HST image, 
as has been employed by Kraj novic et al.1 (12009ft . We used the 
MGE model generated from the WFPC2 79 1W image (dis- 
cussed in which already has the HST PSF removed, and 
convolved it with a model for the OSIRIS PSF. The OSIRIS 
PSF is taken to be three concentric circular 2D Gaussians 
parameterized by the dispersions and relative weights of the 
components. The convolved image is then compared to the 
OSIRIS data cube, summed along the wavelength axis, and 
the parameters of the PSF are varied to determine the best fit. 
We found that the three components have dispersions of 0."01, 
0"04, and 0."28 with relative weights of 0.19, 0.60, and 0.21. 
Thus, we were able to achieve a good AO correction and the 



PSF core has a FWHM of -0." 1 . 

To model the STIS PSF, we used Tiny Tim dKrist & HooU 
120041) to reconstruct the PSF for a monochromatic filter pass- 
band at 865 1 A. Fitting two concentric circular 2D Gaussians 
to the Tiny Tim image, we found that the PSF can be charac- 
terized by a narrow component with a = 0."03 and a weight of 
0.76, plus a broad component with a = 0."12 and a weight of 
0.24. These dispersion values and weights were used as input 
into the stellar dynamical models, which require a PSF to be 
described with circularly symmetric Gaussians. 

For the LRIS PSF, we examined the spatial distribution of 
the eight velocity template stars that were observed on the 
same night as NGC 3998. We fit a single Gaussian to a spatial 
cut at 8600 A through the geometrically rectified 2D spectra, 
and found an average dispersion of a = 0."35. 

7. ORBIT-BASED STELLAR DYNAMICAL MODELING 

Orbit-based dynamic al models using the Schwarzschild 
superposition method dSchwarzschildl [19791) are the stan- 
dard means by which to measure black hole masses from 
stellar kinematics (e.g.. iGebhardt et al.l 12003 1 ; IValluri et alj 
l200l iGuitekin et all l2009at ICappellari et al.l 120091) . They 
are additionall y useful for studying t he orbital structure of 
galaxies (e.g.. |Verolme et al.l l2002bt iRrainovic et al.l 120051 : 
Cappe llari et al.l 120061: I Shen & Gebhardtl 120101). and con 



strain ing dark matter halo propertie s (e.g., [Th omas et al. 



12001 2007. IWeiimans et alj 120091: lForesteir& Gebhardt 



12.010b iMurphy et al.l 1201 ll) . The main strength of this tech- 
nique is that a self-consistent distribution function can be con- 
structed from the observables, without the need for assump- 
tions about the orbital anisotropy. Here, we calculate three- 
integral, triaxial, Schwarzschi l d mo dels using the implemen- 
tation of Ivan den Bosch et alj (120081). w h ich was based upon 
the pre v ious work of iRix et al.l (119971). Ivan der Marel et al.l 
( 1998 al). ICretton et al I (fl999l) . IVerolme & de Zeeuwl (12002a). 
and|Cappellari et al. (2006). We provid e a brief summary of 
the me thod here, but refer the reader to van den Bosc h et al.l 
(2008) for a comprehensive description of model. 

In the model, the galaxy potential consists of contributions 
from the black hole, the stars, and dark matter. We take the 
three-dimensional (3D) stellar mass to be described as the 
sum of multiple coaxial Gaussian distributions, which are de- 
termined by deprojecting the MGE model and assuming a 
constant stellar mass-to-light ratio (M/L). The deprojection is 
carried out by choosing a direction from which the galaxy is 
observed, characterized by the Euler angles 8 (corresponding 
to the inclination angle), tp, and ip. With a set of these view- 
ing angles, each component /' of the MGE model is uniquely 
deprojected into a 3D Gaussian with a shape given by pj, qj, 
and uj, which are the ratios between the intermediate and long 
axes, the short and long axes, and the length of the longest axis 
of the projected Gaussian on the sky to the long axis, respec- 
tively. For details about the deprojection and the relations be- 
tween the viewing angles a nd intrinsic shape parameters see 
Ivan den Bosch et al] (120081) . The Schwarzschild method ad- 
ditionally allows for dark halo properties to be determined, 
but our stellar kinematics extend out to a radius of ~15" (~1 
kpc), and we are unable to differentiate between various dark 
halo models. Consequently, the model results presented in §[8] 
do not include dark matter, and the galaxy potential is calcu- 
lated from the contributions of the black hole and the stars 
alone. However, we do test the effect on the inferred black 
hole mass by incorporating a couple of fixed dark halo mod- 
els in 
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FIG. 6. — A comparison of the stellar kinematics for NGC 3998 measured from the OSIRIS (black circles), STIS (red triangles), and LRIS (blue squares) data. 
The radial velocity, V, the velocity dispersion, a, and the higher-order Gauss-Hermite moments, I13 and /14, which quantify asymmetric and symmetric deviations 
from a Gaussian, are plotted as a function of projected distance from the nucleus. There are two velocity branches in the top left plot due to the galaxy's rotation: 
one side of the galaxy is blueshifted relative to the systemic velocity and the other side is redshifted. The high angular resolution kinematics from both OSIRIS 
and STIS are in agreement with the large-scale data from LRIS . 



Given the galaxy potential, a representative orbit library is 
generated, and the initial conditions of the orbits are chosen 
such that the three integrals of motion are sampled. The orbits 
are integrated and projected onto the observable quantities 
(the location on the sky and the LOS VD parameters) while ac- 
counting for the PSF and aperture binning. Next, weights for 
each orbit are found such that superposition of orbits matches 
the kinematics in the least-squares sense while also satisfying 
contraints imposed by the stellar density in each aperture and 
3D bin. Many dynamical models are calculated with differ- 
ent values for the parameters of interest, namely Mbh, M/L, 
and the viewing angles 9, (j>, and ip (or equivalently the intrin- 
sic shape parameters p, q, and u). Finally, the best model is 
taken to be the one that most closely matches the data in the 
X 2 sense. 

We fit the Schwarzschild models to the OSIRIS, STIS, and 
LRIS kinematics presented in ^5] and use the stellar density 
distribution constructed from the MGE models both with and 
without the innermost Gaussian component discussed in $3] 
The MGE model that includes the central Gaussian compo- 
nent will be discussed in 5(8] and the model that assigns all 
the light from the innermost component to the AGN will be 
described in jj8.ll With four Gauss-Hermite moments mea- 
sured in each of the 90 OSIRIS bins, 8 STIS bins, and 170 
LRIS bins, there are 1072 observables. The orbital libraries 
were set up to sample 25 equipotential shells at radii loga- 
rithmically spaced from 0."006 to 294", with 8 angular and 
8 radial values at each energy to cover the three integrals of 
motion. We bundle together 5 3 orbits with adjacent starting 
positions and sum their observables, resulting in a total of 
600,000 orbits. Such orbital dithering is commonly used in 
the construction of Schwarzschild models to ensure a smooth 
distribution function, and further smoothing can be applied 
after the linear orbital superposition through the adoption of 
a regularization term. Our results presented in §|8]are based 
upon models run without regularization, however we note that 
including a small amount of regularization did not change the 



best-fit black hole mass. 

We estimate the model fitting uncertainty on Mbh (M/L) 
by marginalizin g over M/L (Mbh) an d the shape parame- 
ters. F ollowing Cappella ri et al.l (12009) and iKrainovic et alj 

(2009) , who advise the use of 3cr errors for Mbh measure- 
ments due to the numerical uncertainties associated with 
Schwarzschild modeling, we quote the 3er uncertainties on 
Mbh and M/L for one degree of freedom. This corresponds 
to a change of 9.0 in x 2 from the minimum value. The sta- 
tistical uncertainties on the intrinsic shape parameters are de- 
rived differently, using the c onfidence levels established by 
Ivan den Bosch & van de Venl (120091) and applied during the 
studies of M32 and NGC 3379 bv lvan den Bosch & de Zeeuwl 

(2010) . The confidence interval for the intrinsic shape pa- 
rameters are set based upon the expected standard deviation 
in x 2 , or y/2N b s , where N b s is the number of observables 
used to constrain the model (1072 for NGC 3998). While the 
Mbh determination is influenced by a the innermost kinemat- 
ical measurements, the intrinsic shape parameters are set by 
a much larger number of observables, such that the expected 
scatter in x 2 becomes important and is much larger than the 
standard A% 2 criterion (Ivan den Bosch & van de Ven 2009). 

8. RESULTS 

It is computationally prohibitive to explore parameter space 
for the full range of Mbh, M/L, p, q, and u values simul- 
taneously, so we initially fixed M BH and varied M/L and the 
shape parameters. We sampled axis ratios of 0.60 < p < 1 .00, 
0.40 <q< 0.88, and all possible values of u in steps of 0.06, 
and evaluated 11 /-band M/L values between 4.16 and 6.24 
Mq/Lq. The minimum values of p and q were chosen to 
reflect the smallest values that have been observed in other 
galaxies as well as to probe very triaxial shapes, and the up- 
per bound on q was set by the average flattening of 2D Gaus- 
sians in the MGE model for NGC 3998. With our sampling, 
we have covered 128 different galaxy shapes, 5 of which are 
oblate axisymmetric. This procedure assumes that the shape 
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parameters do not depend on the black hole mass, which was 
shown to be true for the specific ca ses of M32 and NGC 3379 
( Ivan den Bosch & de Zeeuwl 120101) . We tested whether this 
also holds for NGC 3998 by fi rst setting M BH = 2.2 x 10 8 M© 
fthe lde Francesco et al.l ((2006) gas dynamical black hole mass 
measurement adjusted for our assumed distance], and then 
using a larger mass of M B h = 7.9 x 10 8 M© [the mass pre- 
dicted from the Mbh _ relationship when a bulge stellar 
velocit y dispersion of 305 km s" 1 is assumed (Giiltekin et al. 
2009b)]. We found that the best-fit intrinsic shapes differ for 
the two Mbh values tested: an oblate axisymmetric shape, 
with p = 1.00 and q = 0.81 at an effective radius (R e ) is pre- 
ferred when M B h = 2.2 x 10 8 M©, and a round, but triax- 
ial shape, with p = 0.91 and q = 0.81 at 1 R e is found when 
M BH = 7.9 x 10 s M Q . 

Given that the best-fit intrinsic shape changes with Mbh, we 
constructed a model grid by varying M B h between 5.0 x 10 7 
and 5.0 x 10 9 M© and examining 1 1 /-band M/L values be- 
tween 4.2 and 6.2 M©/L©, while also sampling ten shapes. 
These shapes have the lowest ten \ 2 values from the two grid 
runs described above, and they fully encompass the 3a un- 
certainties of p, q, and u when Mbh is set to 2.2 x 10 8 and 
7.9 x 10 s M©. We therefore should be covering the range of 
possible shapes for NGC 3998 during our search for the best- 
fit Mbh and M/L parameters. 

The results of our dynamical models are summarized by 
Figure [7] which displays the contours of % 2 as a function of 
Mbh and M /L after marginalizing over galaxy shape. The best 
model has a \ 2 of 1600, corresponding to a x 2 P er degree 
of freedom (xl) of 1.5, with M BH = (8.1^) x 10 8 M©, I- 
band M/L = 5.0^[j | M© /L©, and an intrinsic shape described 
by p = 0.96 and q = 0.81 at 1 R e . The M BH and M/L er- 
rors represent the 3er statistical uncertainties, and the errors 
for the shape parameters will be discussed below. The kine- 
matics predicted from such a model are compared to the ob- 
served OSIRIS, STIS, and LRIS data in Figures [8] and 
PTOl respectively. We additionally compare the OSIRIS ve- 
locity dispersions to those predicted from models with the 
best-fit mass of M B h = 8.1 x 10 8 M©, a smaller black hole 
with M B h = 1.3 x 10 8 M©, and a larger black hole with 
M BH = 3.6 x 10 9 M© in Fi gure QT| This figure demonstrates 
that our best-fit black hole mass is a reasonable one, and clear 
differences between the predicted velocity dispersions for the 
models with a smaller and larger black hole and the observed 
kinematics can be seen by eye. The more massive black hole 
produces velocity dispersions that are much too large at the 
center, while the smaller black hole is unable to match the 
sharp rise in the observed nuclear velocity dispersion. 

We calculated another model grid holding Mbh fixed at 
Mbh = 8 . 1 x 1 8 M© and allowing M/L and the shape parame- 
ters to vary. This grid allowed us to derive the uncertainties on 
the shape parameters, having previously determined the best- 
fit black hole mass. We again sampled 128 galaxy shapes with 
axis ratios of 0.60 < p < 1 .00, 0.40 < q < 0.88, and all possi- 
ble values of u, and examined 1 1 /-band M/L values between 
4.2 and 6.2 M©/L©. We found that the galaxy shape can be 
described with the ratios p = 0.96^ and q = 0.8 1^ (3a 
uncertainties) at 1 R e . In Figure [T2J we plot the best-fit axis 
ratios p and q and their uncertainties at all radii, extending out 
to 100". We further show the radial variation of the triaxiality 
parameter, T = (1 -p 2 )/(l-q 2 ), on the same plot. An upper 
error bar of zero is measured for q because this ratio is lim- 
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FIG. 7. — Results of the stellar dynamical models after marginalizing over 
the intrinsic shape of the galaxy. At each grey cross, a dynamical model was 
calculated for the specified combination of black hole mass and /-band M/L, 
and the red cross marks the best-fit model. Overplotted are contours of Ax 2 , 
with the inner two contours denoting the Icr and 2<r confidence levels, and 
the thick red contour signifying the 3<r interval for one degree of freedom. 
Contours beyond the 3cr confidence level are separated by a factor of two. 

ited by the observed average flattening of the 2D Gaussians 
from the MGE model. Thus, the best-fit intrinsic shape is as 
round as the HST/CFHT images allow, and is consistent with 
an oblate axisymmetric spheroid. Furthermore, we were un- 
able to place strong constrains on the viewing angles, finding 
that the inclination ranges from 9 = 38° to 9 = 90° (edge-on), 
cf> = -90^ ( 82 , and ip = 90^J (3a uncertainties). These results 
are consistent with previous stellar dynamical studies of other 
early-t ype galaxies, in particular the work of | Krajnovic et alJ 
(2005) and lvan den Bosch & van de Venl (120091) who find that 
the viewing angles are highly degenerate. 

8.1. Additional Sources of Uncertainty 

The 3a errors on Mbh presented above represent the for- 
mal model fitting uncertainty, and account for the random 
noise within the stellar dynamical models. We further ex- 
plored other sources of uncertainty that are not included in 
the statistical errors but that could have an effect the Mbh de- 
termination. We summarize the results of these model grids 
below. 

Stellar Density Distribution: Care must be taken to separate 
the AGN light from the stellar contribution when constructing 
the luminous mass model. Often, this is accomplished by re- 
moving the innermost Gaussian component from the MGE 
model, however NGC 3998 exhibits a steep surface bright- 
ness profile and some starlight may still be contained within 
the central component. In the analysis presented above, we 
assumed that all of the light from the central MGE Gaussian 
was due to the stars. Here, we consider the other extreme, 
where all the light from the inner MGE component comes 
from the AGN. 

After removing the central Gaussian component from the 
MGE model, we followed the same procedure outlined in §|8] 
We first calculated stellar dynamical models fixing Mbh to 
2.2 x 10 8 M© and 7.9 x 10 8 M© in order to determine the most 
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FIG. 8. — The stellar kinematics measured for NGC 3998 from the OSIRIS data (top panels) and the predicted values (bottom panels) from the best-fit stellar 
dynamical model with Mbh = 8.1 X 10 s Mg, /-band M/L = 5.0 Mq /£q, and a shape described by p = 0.96 and q = 0.81 at 1 R e . The velocity map shows that 
the galaxy is rotating rapidly, the velocity dispersion map displays a sharp peak within the inner 0." 1 , and the hi map is anti-correlated with the velocity map. 
The long axis of the IFU was aligned with the major axis of the nuclear gas disk at a position angle of 308° , thus the top of the maps correspond to the northwest 
side of the galaxy. The data and model maps are plotted on the same scale, with the ranges given by the color bar to the right and the minimum and maximum 
values printed at the top of the maps. The kinematics measured from the bins in dark grey were deemed unreliable, and were excluded from the subsequent stellar 
dynamical modeling. 



probable galaxy shapes. Then, we ran model grids sampling 
over the ten best galaxy shapes, while varying Mbh between 
5.0 x 10 7 and 5.0 x 10 9 M© and M/L between 4.2 and 6.2 
Mq/Lq. We found a very round, but oblate, intrinsic shape 
for the galaxy, withM BH = (lO.l^'g) x 10 8 M© and an/-band 
M/L = 4.8+o l M e/ L e ( 3cr uncertainties). This best-fit stel- 
lar dynamical model is a worse description of the observa- 
tions (x 2 = 1619) compared to the best model presented pre- 
viously in ^8](x 2 = 1600). Moreover, when excluding the cen- 
tral MGE component, the dynamical models for every com- 
bination of Mbh and M/L have a higher \ 2 than the stellar 
dynamical models constructed from a mass model with all 12 
Gaussian components. For this reason, and because the black 
hole masses determined using the two MGE models are fully 
consistent within the statistical uncertainties, covering nearly 
identical Mbh ranges, we view the best-fit model from §|8] as 
our final result. 

Dark Halo: Work bv lGebhardt & Thomasl (12009b has raised 
concerns that stellar dynamical black hole mass measure- 
ments may be underestimated if both the black hole and the 
dark matter halo are not simultaneously modeled. The reason- 
ing is that without a dark halo, the M/L may be overestimated 



in order to compensate for the missing mass at large radii. 
The M/L is assumed to be constant in the models, and thus a 
smaller black hole would be needed to match the observed 
nuclear kinematics. By including the contribution of dark 
matter in their re-examination of M87, Gebhardt & Thomas 
(2009) found that the black hole mass increased by about 
a factor of two. Likewise, McConnell et aj] (1201 lal) noted 
a strong dependence of the black hole mass on the dark 
halo when studying NGC 6086. Both objects exhibit shal- 
low stellar luminosity profiles and have massive, concentrated 
dark matter halos, limiting the radial range over which the 
stars dominate the gravitational potential and the assump- 
tion of a constant M/L holds. Moreover, these studies uti- 
lized data in which the quality and/or spatial resolution was 
unable to limit the degeneracy between the black hole and 
stellar mass-to-light ratio in the central regions. In contrast, 
black hole mass measurements based upon high S/N observa- 
tions well within the region influenced by the black hole ap- 
pear to be less sensitive to the inclusion of dark halos during 
the modeling (IShen & Gebhardt! 120 ltit ISchulze & Gebhardll 
l2WltlGebhardt et al.11201 llT " 

We are unable to directly fit for the parameters of a dark 
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FIG. 9. — The observed STIS kinematics (black) and the predicted kine- 
matics (red line) from the best-fit dynamical model. The velocity, velocity 
dispersion, /13, and h$ are plotted as a function of position, relative to the 
nucleus, along the STIS slit. The slit was aligned with the galaxy's major 
axis. We were unable to measure the stellar kinematics within ~0''2 of the 
nucleus due to AGN contamination, and model predictions were not gener- 
ated for this region. The kinematics are similar to those measured from the 
OSIRIS data. 



halo model because our kinematics do not extend out to 
large enough radii. Instead, we selected a fixed model 
for the dark halo. Two commonly used dar k halo models 
are th e Navarro-Frenk- White (NFW) profile (tNavarro et al.l 
11996 ] ) and a distribution based on a cored logarithm ic poten- 



tial (iBinney & Tremainelll987t iThomas et al J 120051) . Previ- 
ous work ha s found that both parameterizations give consis- 
tent results (IThomas et al.1 120071: iGebhardt & Thomasl 120091: 
iMcConnell et alj|201 lal) . and here we adopt the profile from 
a logarithmic potential given by 



3r 2 r +r 2 



4vrG Of + r 2 ) 



2\2 ' 



(1) 



The parameters V c and r c are the asymptotically constant cir- 
cular velocity and the core radius, within which the dark 
matter density is approximately constant. We fixed V c and 
r c to values of 407 km s" 1 and 10.7 kpc, respectively, 
which were chosen for NGC 3998 using the galaxy B- 
band luminosit y of 1.6 x 10 10 L reported by the Hyper- 
led a da tabase ( Patur elet al.ll2003l) and the empiri cal relations 
by IThomas et alJ ( 120091) . The relationships in Thomas et al. 
(2009) provide a way to select physically motivated param- 
eters for the dark matter h alo, and was similarly used by 
ISchulz e & Gebhardt d201 ll) to test the effects of dark matter 
on the black hole mass. 

When including the dark halo in our stellar dynamical mod- 
els, we did not observe a significant change in the black hole 
mass, measuring Mbh = x 10 8 M and an /-band 

M/L = 5.0^0 4 M Q /L Q (3cr errors). The decrease in the best- 



fit black hole mass is in the opposite direction as anticipated, 
but is well within the statistical uncertainties and is likely due 
to noise in the dynamical models. A similar result was seen 
when we incorporated a dark halo that was twice as massive 
within the radial range covered by the large-scale LRIS kine- 
matics. Therefore, neglecting the dark matter contribution 
during the dynamical modeling presented in ^8] did not lead 
to a biased Mbh determination. We attribute this result to the 
high quality OSIRIS data, which probes regions well within 
the influence of the black hole and reduces the degeneracy 
between Mbh and M/L, as well as to the galaxy's strong rota- 
tion, which also makes the measurement less sensitive to the 
dark halo. 

Central OSIRIS Kinematics: Contributions from a non- 
stellar nucleus dilute absorption lines and make the mea- 
surement of reliable stellar kinematics challenging. In some 
instances, the AGN contribution is so strong that the stel- 
lar features are no lo nger visible in the nuclear spectra 
dCappellari et al.ll2009t) . This is not the case for NGC 3998, 
however we do measure large errors for the kinematics ex- 
tracted from four of innermost OSIRIS lenslets. As an exam- 
ple, at these locations we found velocity dispersions between 
500 and 560 km s" 1 with error bars of 42 - 146 km s -1 . Be- 
yond these four lenslets, but still within ~0."1 of the nucleus, 
we continued to measure large velocity dispersions of a ~ 450 
km s" 1 but with more reasonable errors of ^20 km s -1 . The 
large uncertainties from the inner OSIRIS bins indicate that 
the pPXF routine could not model the AGN contribution in 
a consistent fashion within 0."05 from the nucleus, which 
in turn may have biased these kinematical measurements. 
Therefore, we excluded the measurements from the central 
five OSIRIS bins when fitting the orbit-based stellar dynam- 
ical models. Doing so led to a best-fit model with Mbh = 
(8.1 ±2.0) x 10 s M and /-band M/L = 5.0^ M /L . The 
values and 3er uncertainties for Mbh and M/L are the same 
as found previously in $8] when all of the OSIRIS kinematics 
were fit during the modeling. Even if the kinematics measured 
from four of the innermost OSIRIS lenslets are biased due to 
AGN contamination, there appears to be no subsequent effect 
on the black hole mass. This result is not entirely surprising 
because the stellar dynamical models take into account the er- 
rors on the input kinematics, and we have a number of other 
reliable measurements of the LOSVD within the black hole 
sphere of influence. 

OSIRIS PSF: As discussed in igj] the OSIRIS PSF was de- 
termined by comparing the deconvolved HST WFPC2 F79 1 W 
image to the OSIRIS /T-band data cube summed along the 
wavelength axis. This approach assumes that the comparison 
image has a higher resolution than the collapsed data cube 
and that there is no strong color gradient. Given the difficul- 
ties associated with measuring the PSF, we therefore addition- 
ally considered two extreme parameterizations of the OSIRIS 
PSF to quantify changes in Mbh- Both PSFs were composed 
of two concentric, circular, 2D Gaussians, and the dispersion 
of the second component was set at 0"28. We chose 0"28 
because this is the dispersion of the broad component in the 
OSIRIS PSF model presented in ^6] We then constructed a 
"poor" PSF, such that the narrow component had a dispersion 
of 0"084 and contributed 15% of the flux to the total PSF, and 
a "good" PSF, whose narrow component had a dispersion of 
0."025 (roughly the ZT-band diffraction limit for a 10m tele- 
scope) and contributed 95% of the flux to the total PSF. The 
stellar dynamical models based upon the "bad" PSF yielded 
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FIG. 10. — The observed LRIS kinematics (black) and predicted values from the best-fit dynamical model (red line). The LRIS observations were obtained at 
four position angles: along the major axis of the nuclear gas disk (PA = 308°), along the minor axis (PA = 218°), and at two intermediate angles (PA = 353° and 
PA = 263°). For each of the slit positions, the velocity, velocity dispersion, /13, and /14 are plotted as a function of location along the slit relative to the nucleus. 
These large-scale kinematics exhibit similar features to those seen from the high angular resolution data, which include rapid rotation, a steep rise in the velocity 
dispersion, and an anti-correlation between /13 and V. Furthermore, A4 is symmetric about the center with a slight peak at the nucleus. 
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FIG. 1 1. — Comparison between the observed velocity dispersions over the OSIRIS field of view and the predicted ones for various dynamical models. The 
models were constructed with different black hole masses (the best-fit Mbh of 8.1 X 10 8 Mq, a larger black hole with Mbh = 3.6 X 10 9 Mq, and a smaller black 
hole of Mbh = 1.3 X 10 Mq), but have the same /-band M/L (5.0 Mq/Lq) and intrinsic shape (p = 0.96 and q = 0.81 at 1 R e ). All maps are plotted on the same 
scale, given by the color bar to the right along with the minimum and maximum values listed on top of the left-most panel. Observed kinematics measured from 
the grey bins were deemed unreliable, and we do not show the model predictions in these bins. 
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FIG. 12. — The intrinsic shape of NGC 3998. The ratio, p, between the 
intermediate and long axes (solid black line), the ratio, q, between the short 
and long axes (dash-dotted red line), and the triaxiality parameter T (dashed 
blue line) are plotted as a function of distance from the nucleus. The 3er un- 
certainties are given by the range covered by the black backward, red vertical, 
and blue forward slashes for p, q, and T, respectively. The kinematical mea- 
surements from the large-scale LRIS data extend out to a radius of ~15". 

a best-fit model with M BH = (8.7t^) x 10 8 M and /-band 
M/L = 5.0!^ Mq/Lq, and the models using the "good" PSF 
produced best-fit parameters of M BH = (8.I+I4) x 10 8 M and 
/-band M/L = 5.0 ± 0.2 Mq/Lq. Thus, even for two extreme 
PSF models, the black hole mass remains within the 3a sta- 
tistical uncertainties presented in §|8] 

Symmetrizing the Kinematics: We removed the systematic 
offsets in the odd Gauss-Hermite moments (e.g., the galaxy's 
systemic velocity) with th e poin t-symmetrization code of 
Ivan den Bosch & de Zeeuwl d2010). In addition to this step, 
the common practice is to make further adjustments to the 
kinematical measurements and their errors so that the veloc- 
ity fields become bi- or point-symmetric, thereby reducing 
the observational noise. However, the symmetrization pro- 
cess assumes that the galaxy nucleus is centered on one of 
the lenslets, which is not necessarily true for the OSIRIS 
data. Therefore, in §|8] we calculated stellar dynamical mod- 
els without first requiring the observed velocity fileds to be 
symmetric. Instead, if we do symmetrize the input kinemat- 
ics, we find a best-fit model with M BH = (9.5if 2) x 10 8 M Q 
and an /-band M/L = 4.8^2 Mq/Lq. As expected, the best 
model is a much better description of the (symmetrized) ob- 
served kinematics, with \ 2 = 459 and xl = 0.43, and larger 
statistical 3a uncertainties on M B h are found. Although the 
best-fit black hole mass increased by 17% from the mass de- 
rived previously, it remains within the model fitting uncertain- 
ties determined in §|8] and we continue to prefer the dynamical 
model constrained with the non-symmetrized kinematics. 

9. DISCUSSION 

Through orbit-based stellar dynamical modeling we mea- 
sured a mass of M BH = (8.1!^) x 10 s M Q for the black hole 
in NGC 3998. We report 3a statistical errors on the black hole 
mass, but ran additional model grids to assess the robustness 
of our mass measurement. We found no significant changes 
(outside the random noise of the models) to the black hole 



mass due to these sources of uncertainty, and so we do not 
make any additional adjustments to the black hole mass or 
its error bars. With this black hole mass, and a bulge stellar 
velocity dispersion of 272 km s" 1 (see j39.3| below), the black 
hole sphere of influence is r sp h e re = 0"7, which is well resolved 
by the observations. We now examine the orbital structure of 
galaxy, compare the stellar dynamical black hole mass mea- 
surement to the existing gas dynamical mass, and place NGC 
3998 on the M B h - host galaxy relationships. 

9.1. Orbital Structure 

We examined the internal orbital structure of NGC 3998 
using the orbital weights found with the Schwarzschild su- 
perposition method for our best fitting model (presented in 
®. Defining the tangential velocity dispersion as a} = 
(o^ + Og)/2, with (r,8,(f>) being the usual spherical coordi- 
nates, we plot the ratio a,-/ a, as a function of radius in 
the top panel of Figure [13] Near the nucleus, the velocity 
dispersion becomes radially anisotropic reaching values of 
a r /a, 1.5, while the velocity dispersion is isotropic, de- 
viating by at most ~10% from a r /a, = 1, at radii beyond 
^0"1. The radial anisotropy near the galaxy's center can 
be attributed to the large fraction of box orbits: we found 
that the box orbits contribute 60 - 80% to the stellar mass 
within the inner ~0"1 (see bottom panel of Figure [L31. Sim- 
ilarly large contributions from box orbits have been seen 
when appl ying triaxial Schwarzschild models to other objects 
as w ell dWeiimans et alj 120091 : Ivan den Bosch & de Zeeuw 
2010). Outside of r sp h e re> our best-fit model is made up of 
^65% short-axis tube orbits and ^20% long-axis tube orbits. 

NGC 3998 is classified as an SO galaxy, and indeed, our best 
fitting dynamical model shows evidence for both a bulge and a 
disk component. In Figure [14] we show the mass distribution 
along orbits as a function of average radius and spin, \ z . The 
spin is defined as \ z = J z x (r/ a), where J z is the average an- 
gular momentum along the z-direction, f is the average radius, 
and a is the average second moment of the orbit. The figure 
shows the clear presence of a non-rotating bulge component 
(-0.2 < A. < 0.2) and a rotating disk component (\ z > 0.2), 
making up 52% and 39% of the mass within the radial range 
covered by our kinematic measurements. The remaining 9% 
of the mass comes from orbits with A z < -0.2 that fall outside 
the radial range covered by the OSIRIS kinematics. Thus, the 
near isotropy at radii larger than ~0" 1 seen in FigureQ~3]is due 
to the presence of a bulge component, while the disk causes 
the strong rotation that is seen in the OSIRIS, STIS, and LRIS 
kinematics in Figures [8] [9] andflOl 

9.2. Comparison to the Gas Dynamical Measurement 

By modeling the HST STIS gas kinematics as a thin, circu- 
larly rotating disk, Ide Francesco et all (|2006) derived a black 
hole mass of M BH = (2.7^'q) x 10 8 M Q (2a uncertainties). 
Scaling the black hole mass for the distance adopted here, 
their mass measurement becomes M B h = (2.2+[g) x 10 8 M Q . 
The stellar dynamical black hole mass is therefore signifi- 
cantly larger, by nearly a factor of four, compared the gas 
dynamical measurement. We cannot resolve the discrepancy 
between the stellar and gas dynamical mass measurements for 
the black hole in NGC 3998, but discuss a few of the uncer- 
tainties associated with both methods below. 

The gas kinematics measured from the STIS data clearly 
show organized rotation. However, the line widths pre- 
dicted from a rotating, dynamically cold, thin-disk model 
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FIG. 13. — Orbital structure of NGC 3998 from our best fitting stellar dy- 
namical model. The anisotropy, a,/ a,, (top) and the orbit type (bottom) are 
plotted a s function of radius, covering the extent of the kinematic measure- 
ments. The vertical grey dot-dashed line designates the black hole sphere of 
influence. NGC 3998 is mostly isotropic (indicated by the horizontal dashed 
grey line in the top panel) but shows a radial bias within ~0" 1 . Near the 
nucleus, box orbits (black dotted line) dominate the galaxy, while short-axis 
tube orbits (red solid line) become important at larger radii. Long-axis tube 
orbits are shown by the blue dashed line. 
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FIG. 14. — Distribution of mass along the orbits from our best fitting dy- 
namical model as a function of average radius and spin. The radial range 
extends over the region covered by our kinematic measurements. A non- 
rotating bulge (-0.2 < A- < 0.2) and rotating disk (A- > 0.2) component are 
clearly present. 

do not fully ma t ch th e observed velocity dispersions. 
Ide Francesco et all (12006b found that the model underesti- 
mates the observed rise in the nuclear velocity dispersion by 
^20%, although larger inconsistencies were seen in the two 
slit positions immediately adjacent to the central slit. With- 
out re-analyzing the STIS data for NGC 3998, it is difficult 
to determine whether the addition of an intrinsic velocity dis- 
persion to the models would greatly improve the agreement 
between the predicted and observed line widths, and by how 



much the gas dynamical black hole mass measurement would 
then increase if the excess velocity dispersion was considered 
dynamically important. Previous gas dynamical studies that 
have assigned a dynamical origin to the internal velocity dis- 
persion have found varying degrees to which the black hole 
mass is affected. For example, using an asymmetric drift cor- 
rection to estimate the circular velocity f rom the observed 
mean rotational speed, iBarth et ail (120011) found a 12% in- 
crease in the best- fit black hole mass for NGC 3245 while 
Wals h et alj (|2010) found that the mass of the black hole in 
M84 doubles. However, we note that the differences between 
the observed line widths and the predictions from a dynami- 
cally cold, rotating disk model for NGC 3245 and M84 were 
larger than those observed for NGC 3998. Admittedly, the 
factor of four increase needed to make the best-fit gas dynam- 
ical mass to match the best-fit stellar dynamical measurement 
is quite large, and the discrepancy probably cannot be entirely 
resolved by assigning a dynamical origin to an intrinsic veloc- 
ity dispersion. 

Another source of uncertainty in the gas dynamical black 
hole s mass comes from the in clination of the emission-line 
disk. Ide Francesco et al.l (f2006) found that all models with in- 
clinations below 70° produced equally acceptable fits to the 
data at the 2a level, and that the inferred black hole masses 
differed by about an order of magnitude for inclination an- 
gles i= 10° and 70°. A similarly large difference in the NGC 
3998 black hole mass fo r extre me inclination angles was re- 
ported by Beifiori et al. (2009) , who estimated masses of 
M BH = 3.4 x 10 8 M Q and 4.9 x 10 7 M Q ( adjusted to our dis- 
tance) for inclinations of i = 33° and 88°. Ide Francesco et all 
(2006) ruled out very face-on orientations, i < 21°, because 
such angles produced V-band mass-to-light ratios that are too 
high compared the predictions of a Salpeter initial mass func- 
tion (IMF) for an old stellar population. They ultimately re- 
ported results using a disk inclination of 30°, and further com- 
bined their statistical uncertainties on Mbh with the uncer- 
tainties resulting from the allowed range in inclination angles 
between 27° and 70° to derive the 2cr error bars on Mbh. In- 
terestingly, the results of Cappellari et al] (120121) indicate that 
the most massive galaxies have M/L values that are higher 
than the predictions of a Salpeter IMF for an old stellar pop- 
ulation, thus lower inclinations of the gas disk in NGC 3998 
may be acceptable, which would then lead to an increase in 
the gas dynamical mass measurement. 

When constructing the stellar dynamical models presented 
here, we did not assume axisymmetry or a specific viewing 
orientation. However, it is not feasible to explore all possi- 
ble values for the five parameters simultaneously, and we em- 
ployed a method that begins by holding fixed Mbh while vary- 
ing the other parameters to find the best-fit intrinsic galaxy 
shape before searching for the best-fit Mbh- Such an approach 
does not guarantee that a global minimum will be found, and 
we attempted to alleviate this problem by sampling over ten 
galaxy shapes (instead of just the best-fit shape) while vary- 
ing Mbh and M/L. The ten shapes were taken from the ten 
best models when Mbh was held at two different values, and 
the shapes cover the 3cr ranges of p, q, and u for the two 
black hole masses. Thus, we should have covered the pos- 
sible galaxy shapes while searching for the best fitting Mbh 
and M/L parameters. 

One assumption that is made when constructing the stellar 
dynamical models is that the M/L does not vary with radius. 
NGC 3998 has both bulge and disk components, and could 
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also contain a nuclear star cluste r ([Gonzalez Delgado et al.l 
2008; Neuma veF& Walcheril20l2l) . If each of these compo- 
nents has a very different stellar population, then a single M/L 
is not a realistic description of the system. Recent work (e.g., 
Ivan Dokkum & Conroyl l20Toh has advocated for a bottom- 
heavy IMF in massive early-type galaxies. Thus, if the bulge 
component of NGC 3998 is indeed rich in dwarf stars, the 
larger M/L for the bulge region would lead to a smaller stellar 
dynamical black hole mass measurement. However, allowing 
for a radial variation in M/L, in addition to sampling the Mbh 
and three shape parameters, is too computationally expensive 
to be attempted here. In contrast, the gas dynamical black 
hole mass measurement, which relies on the gas kinematics 
at small radii where the stellar potential is dominated by the 
bulge, is insensitive to M/L variations. While our assump- 
tion of a constant M/L, may effect the derivation of the stellar 
dynamical black hole mass, the other systematic effects of- 
ten associated with the stellar dynamical method appear to be 
minimal. For example, we found the black hole mass was in- 
sensitive to the inclusion of a of dark halo in the model, as 
well as to the exact OSIRIS PSF that was adopted. Moreover, 
we found good agreement between the kinematics measured 
from different instruments using different template libraries, 
suggesting that template mismatch is not a dominate source 
of uncertainty. 

9.3. The Black Hole Scaling Relations 

When placing NG C 3998 on the M BH - <r* relation, 
iGiiltekin et al.l (l2009bl) adopted a bulge stellar velocity dis- 
persion of 305 km s . The value was taken from the Hy- 
perleda database, and was based upon the dispersion values 
found in the literature at that time, however we are in a posi- 
tion to directly measure the effective stellar velocity disper- 
sion (er e ) of the bulge component from our data and mod- 
els. Recent studies of very massive black holes have elected 
to exclude data wit hin r S ph Rre when measuring the effective 
velocity dispersion (j Gebhardt et al.l 1201 ll Uardel et al.l 1201 lb 
McConnell et all 1201 lbh . and here we follow the same ap- 
proach. We weight the a and V measurements obtained from 
the LRIS major axis slit with th e surface brigh t ness pr ofile 
set by our MGE model, following Giilte kin et al.l d2009bt) . but 
integrate from r sp h e re to the bulge effective radius. The bulge 
effective radius is rather uncertain for NGC 399 8, with lit- 
erature values ranging bet ween 4. "7 < R e < 18 "3 (ISani et alJ 
1201 U iBaggett et al] [l998). which in turn produces some un- 
certainty in the measurement of cr e . Taking an average of 
the R f measurements determined from recen t optical images 
(lFisheretal.ll 19961: ISanchez-Portal et al.l [20041) . we find that 
R e =10"7 and a e = 272 km s" 1 for NGC 3998. For compari- 
son, effective stellar velocity dispersions of er e = 282 km s" 1 
and 270 km s" 1 are measured for R e =4 ."7 and 18. "3, respec- 
tively. We note, however, that there is just a single measure- 
ment of the large-scale stellar kinematics on each side of the 
nucleus beyond a radius of 12", making it difficult to deter- 
mine cr e for a bulge effective radius of 18. "3 from the LRIS 
data. 

These a e measurements are based upon the large-scale 
kinematics extracted from the major axis LRIS slit, however 
the best-fit stellar dynamical model from fallows us to pre- 
dict the luminosity-weighted second moment from a circular 
aperture of radius R e . We measured an effective velocity dis- 
persion of <7 e = 239 km s" 1 when a bulge effective radius of 
R e =10"7 is used and the central regions within r sp hei- e are ex- 



cluded. This method provides a formally more correct esti- 
mate of cr e , as a circular aperture is being considered instead 
of relying on a single slit position, but most of stellar veloc- 
ity dispersions on the Mbh _ relationship w ere derived us- 
ing lon g-slit data and the definition given by Giiltekin et al. 
(2009b). Thus, for consistency purposes we continue to con- 
sider the measurements of a e made from the LRIS data. 

With an effective stellar velocity dispersion of <r e = 272 
km s" 1 , the black hole mass predicted from the Mbh _ 0* 
relation is 9.4 x 10 8 M p, using the recent calibration of 
iMcCon nell et al. (2011b). Instead, black hole masses of 
6.5 x 10 8 M Q and 4.9 x 10 8 M m are predicted from the 
Mbh — cr* correlation c alibrated by Gra ham et al.l (1201 ll) and 
Giilte kin et al.l (l2009bl) . Thus, our stellar dynamical black 
hole mass measurement falls within the expectations of the 
relationship when calculating a e using the large-scale LRIS 
kinematics and a bulge effective radius of 10. "7. 

To place NGC 3998 on the M B h - Ami correlation, we use 
the total, extinction corrected, V-band luminosity of 9.7 x 
10 9 Lq for NGC 3998 f rom the Third Reference C atalog 
of Bright Galaxies (RC3; de Vaucoule urs et al.l [199 1 ). The 
average of the bul ge-t o-total ratio (B/T) va l ues r eported by 
iFisher et all (119961) and lSanchez-Portal et al l ( 120041) is B/T = 
0.77, which suggests a V-band luminosity of 7.5 x 10 9 L© 
for the bulge of NGC 3998. Such a bulge luminosity trans- 
lates into a prediction of 7.2 x 10 7 M and 5.0 x 10 7 M Q 
for the black hole mass using the V-band M B h - Lp u \ reb - 
tion of IMcConnell et alJ d2011bl) and IGiiltekin et all d2009b). 
Even when deliberately overestimating the bulge luminosity, 
and using the total, V-band luminosity reported in RC3, the 
MfiH-^hni relation sugge sts a black hole mass of 9.6 x 1O 7 M 
(McConnell et al. 2011b). The stellar dynamical measure- 
ment therefore places NGC 3998 well above the Mbh _ ^bui 
relationship. In Figure[l5] we show the location of NGC 3998 
on the black hole mass scaling relations. 

10. CONCLUSIONS 

To conclude, we have studied the stellar dynamics of the 
nearby, SO galaxy, NGC 3998. From Keck LGS AO OSIRIS 
observations in the K band and archival HST STIS data cov- 
ering the Ca II triplet lines, we mapped out the 2D kinematics 
within ^2"of the nucleus with a superb angular resolution of 
~0."1, thereby resolving the gravitational sphere of influence 
of the black hole. In addition, we obtained long-slit data at 
four PAs with Keck LRIS to measure the large-scale stellar 
kinematics, extending to ~1 R e , which is essential for con- 
straining the stellar orbital distribution. We found that the 
galaxy is rapidly rotating with V ~ ±200 km s" 1 and exhibits 
a very sharp rise in the velocity dispersion to values of a ^500 
km s" 1 . Our high-quality spectroscopic data further allowed 
us to quantify the LOSVD's asymmetric and symmetric devi- 
ations from a Gaussian through the /13 and /14 Gauss-Hermite 
moments. We observed an anti-correlation between /13 and V, 
and slight peak in at the center from the large-scale data. 

Combining the kinematics with the luminosity density mea- 
sured from a HST WFPC2 F791W image and a CFHT 
WIRCam K s image, we constructed three-integral, triaxial, 
Schwarzschild models. The intrinsic shape of the galaxy is 
very round (as round as the imaging observations allow), and 
oblate, with axis ratios p = 0.96^ ^ and q = O.8I+033 at a ra- 
dius of 10 "7. Although we are unable to place strong con- 
straints on the shape parameters/viewing orientation, by sam- 
pling over them, we were able to extract a robust measurement 
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FIG. 15. — Locatio n of NGC 3998 on the M BH - a* (left) and V-band M BH -Lbul (right) relationships. The correlations calibrated by Giiltekin et al. (2009b), 
IGraham et ail d20TTb . and McConnell eTafl <201 lbf) are given by the dotted, dashed, and dot-dashed lines, respectively. The errors on the bulge stellar velocity 
dispersion show the range of values that were calculated using the large-scale LRIS data between r sp hc re and the largest and smallest bulge effective radius 
measurements found in the literature [R e =18'' 3 iBaggett et al. 1998) and R e =4"7 iSani et al. 2011 )]. Simila rly, the errors o n the bulge luminosity wer e 
determined by applying the smallest and largest values of B/T found in the literature [B/T = 0.41 ( Sani et al. 20 11) and B/T = 0.83 I Sanchez-Portal et al. 2004)] 
to the total V-band luminosity given in RC3. 



of the black hole mass and /-band stellar mass-to-light ratio, 
finding M BH = (8.1!^) x 10 8 M Q and M/L = 5.0^ M /Lq. 
Additional model grids were run to assess possible systematic 
effects. We tested the effect on the black hole mass when as- 
signing all the light from the central MGE component to the 
AGN, incorporating a fixed dark halo model, excluding the 
central OSIRIS kinematics that may be biased due to AGN 
contamination, adopting two extreme models for the OSIRIS 
PSF, and symmetrizing the kinematic measurements. We did 
not see any significant changes (outside of the modeling fit- 
ting uncertainties) to the best-fit Mbh or M/L. 

With the stellar dynamical black hole mass measurement, 
NGC 3998 is consistent with Mbh _ o* when using a bulge 
stellar velocity dispersion of 272 km s" 1 , but well off the 
MBH-^bui correlation, even when overestimating the bulge 
luminosity. Also, the stellar dynamical black hole mass mea- 
surement is larger than the existing gas dynamical measure- 
ment, with the masses differing by close to a factor of four. 
NGC 3998 is now one of eight galaxies for which both stel- 
lar and gas dynamical modeling have been used to measure 
the mass of the central black hole. However, the gas kinemat- 
ics turned out to be strongly disturbed in two of the galaxies 
and a clear best-fit mass and the associated uncertainties could 
not be measured for another object, preventing a meaningful 
comparison between the stellar and gas dynamical techniques. 
The five remaining comparison studies have produced mixed 
results, ranging from excellent agreement between the two 
mass measurements to the stellar dynamical black hole mass 
exceeding the gas dynamical determination by a factor of five. 
Of the three cases in which the black hole masses disagree, the 
stellar dynamical measurement is always larger than the black 
hole mass derived with gas dynamical modeling. With such 
a limited sample, clearly more cross-checks are necessary be- 
fore conclusions can be made regarding the consistency of the 
gas and stellar dynamical techniques, the subsequent effects 
on the Mbh _ host galaxy relationship, and the magnitude and 
distribution of the cosmic scatter in the correlations. Target- 
ing galaxies with dust lanes could be particularly useful, as 



well-ordered, symmetric dust l anes suggest the presence of a 
regularly rotating gas disk dHo et al.l 12002b and can also be 
used to constrain galaxy inclinations. 
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